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ABSTRACT 

We show that the discrepancy between the TuUy-Fisher relation and the luminosity 
function predicted by most phenomenological galaxy formation models is mainly due 
to overmerging of galaxy haloes. We have circumvented this overmerging problem, 
which is inherent in both the Press-Schechter formalism and dissipationless N-body 
simulations, by including a specific galaxy halo formation recipe into an otherwise 
standard N-body code. This numerical technique provides the merger trees which, to- 
gether with simplified gas dynamics and star formation physics, constitute our imple- 
mentation of a phenomenological galaxy formation model. Resolving the overmerging 
problem provides us with the means to match both the I-band Tully-Fisher relation 
and the B and K band luminosity functions within an £7 = 1 sCDM structure forma- 
tion scenario. It also allows us to include models for chemical evolution and starbursts, 
which improves the match to observational data and renders the modelling more real- 
istic. We show that the inclusion of chemical evolution into the modelling requires a 
significant fraction of stars to be formed in short bursts triggered by merging events. 

Key words: cosmology: theory - dark matter - large-scale structure of Universe - 
galaxy formation 



> 1 INTRODUCTION 



There has been significant recent progress in the study of 
galaxy formation within a cosmological context, mainly due 
to a phenomenological or 'semi-analytical' approach to this 
problem. The idea is to start with a structure formation 
model that describes where and when galactic dark haloes 
form. A simple description of gas dynamics and star for- 
mation provides a means to calculate the amount of stars 
forming in these haloes. Stellar population synthesis mod- 
els then provide the spectral evolution, i.e. luminosities and 
colours, of these galaxies. 

Many physical processes are modelled as simple func- 
tions of the circular velocity of the galaxy halo. Therefore, 
the Tully-Fisher relation is the most obvious observational 
relation to try and predict, as it relates the total luminos- 
ity of a galaxy to its halo circular velocity. However, most 
phenomenological galaxy formation models do not simulta- 
neously fit the I-band Tully-Fisher relation and the B or K 
band luminosity function. When one sets the model parame- 
ters such that the Tully-Fisher relation has the right normal- 
ization, the luminosity functions generally overshoot (e.g. 
Kauffman, White & Guiderdoni 1993; Kaufmann, Colberg, 
Diaferio & White 1998), certainly for the = 1, i/o = 50 
km s~^ Mpc~^ standard CDM cosmology (in the form given 



by Davis et al. 1985) that we consider in this paper. Alter- 
natively, when making sure that the luminosity functions 
matches by changing some of the model parameters, the 
Tully-Fisher relation ends up significantly shifted with re- 
spect to the observed relation (e.g. Cole et al. 1994; Heyl et 
al. 1995). 

In order to keep the modelling as analytical as possible, 
an extension to the Press & Schechter (1974) prescription 
for the evolution of galaxy haloes (e.g. Bond et al. 1991; 
Bower 1991; Lacey & Cole 1993; Kauffman & White 1993) 
has been a popular ingredient for implementations of a phe- 
nomenological theory of galaxy formation. The extended 
Press-Schechter (EPS) formalism keeps the galaxy forma- 
tion model as analytical as possible, and allows fast realiza- 
tions of halo populations. However, this prescription has the 
disadvantage that the large-scale phase-space distribution of 
haloes is unknown, as is any information on their internal 
structure and kinematics. For many applications this in- 
formation will be of great use. Furthermore, the properties 
and formation histories of individual haloes do not match 
the corresponding properties and histories as predicted by 
the EPS formalism (e.g. Lacey & Cole 1993; White 1996). 

More importantly, the EPS formalism is designed to 
identify collapsed systems, irrespective of whether these con- 
tain surviving subsystems. This 'overmerging' of subhaloes 
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into larger embedding haloes produces a top-heavy galaxy 
halo mass function, and is therefore relevant to the problem 
of matching both the luminosity function and the Tully- 
Fishcr relation. Traditional N-body simulations suffer from 
a similar overmerging problem (e.g. White 1976), which is 
of a purely numerical nature, caused by two-body heating 
in dense environments (Carlberg 1994; van Kampen 1995) 
when the mass resolution is too low. This is why the statis- 
tical properties of EPS haloes and haloes found in N-body 
simulations still agree quite well (e.g. Efstathiou et al. 1988; 
Somerville et al. 1998 and references therein), despite the 
overmerging. 

In order to circumvent these problems, we use an N- 
body simulation technique that includes a built-in recipe for 
galaxy halo formation, designed to prevent overmerging (van 
Kampen 1995, 1997), to generate the halo population and 
its formation and merger history. We show that this alone 
significantly changes the predictions for the TuUy-Fisher re- 
lation and the luminosity function. In fact, it resolves most 
of the discrepancy sketched above, and allows us to make 
the modelling more realistic by adding chemical evolution 
and a merger-driven bursting mode of star formation to the 
modelling. 

Star formation and its feedback to the interstellar 
medium are the least understood ingredients of theories of 
galaxy formation. Usually one considers a single mode of 
star formation, being cither a continuous but decaying rate 
of star formation from gas that cools within dark haloes, 
or a discontinuous, bursting mode of star formation where 
stars are formed at a high rate in short periods. Here we 
will consider both, where the continuous mode of star for- 
mation encompasses the case of many short bursts of star 
formation. We reserve the term 'starburst' for single major 
starburst events, lasting on the order of 0.1 Gyr, as ob- 
served in luminous infrared galaxies, for example (Sanders 
& Mirabel 1996). 

Once stars arc formed, wc apply the stellar population 
synthesis models of Jimenez et al. (1998) to follow their evo- 
lution. We have enhanced these models with a model for 
the evolution of the average metallicity of the population, 
which depends on the starting metallicity. Feedback to the 
surrounding material means that cooling properties of that 
material will change with time, affecting the star formation 
rate, and thus various other properties of the parent galaxy. 

Several of our ingredients are thus different from earlier 
work; we use numerical simulations for the formation and 
evolution of the galaxy haloes, and we use newer population 
synthesis models, which include chemical evolution. Nev- 
ertheless, it is still useful to start from a model that is as 
similar as possible to a published one. This allows us to as- 
sess the changes in the model results due to the differences 
in modelling. We chose use the model of Cole et al. (1994) 
as a starting point. 

The lay-out of this paper is as follows: we first discuss 
the main problems in Section 2, focusing on the overmerging 
of galaxy haloes, and ways of resolving the problems. In 
Section 3 we generate distributions of galaxies for which we 
calculate luminosity functions and the TuUy-Fisher relation. 
We compare these to those found in previous work, and to 
observed ones. Wc discuss the interpretation of the results, 
pitfalls, and possible extensions, in the final section. 



2 THE PHENOMENOLOGICAL MODEL 
FOR GALAXY FORMATION 

In this section, we summarize how gala^xy formation is mod- 
elled in this work. In many respects, we deliberately follow 
the approach of other groups, so as to aid comparison of 
results. Nevertheless, our methods differ in some important 
respects, as outlined below. A detailed discussion of the 
modelling is given in Appendix A; this section is intended 
as an overview for non-specialists. 

2.1 Overview 

The key ingredients of the model are: 

(1) The merging history of dark-matter haloes. This is often 
treated by Monte-Carlo realizations of the analytic 'ex- 
tended Press-Schechter' formalism. That formalism is 
a successful description of collapsed systems with den- 
sities of order 200 times the mean, and of the history of 
merging that produced such systems. However, the the- 
ory ignores substructure: a cluster-scale halo is treated 
as a single system. The alternative is to measure halo 
substructure directly via N-body simulation. Special 
techniques are used in this paper to prevent galaxy- 
scale haloes undergoing 'overmerging' owing to inade- 
quate numerical resolution. By allowing substructure to 
survive, the merger history is simplified and the number 
of halo mergers is reduced. 

(2) The merging of galaxies within dark-matter haloes. 
Each halo contains a single galaxy at formation. When 
haloes merge, a criterion based on dynamical friction 
is used to decide how many galaxies exist in the newly 
merged halo. The most massive of those galaxies be- 
comes the single central galaxy to which gas can cool, 
while the others become its satellites. This approxi- 
mate argument is one of the more uncertain parts of 
phenomenological galaxy formation schemes. We still 
use the dynamical friction method, but our treatment of 
halo overmerging means that part of the galaxy merg- 
ing process is treated directly. This should make our 
results more robust. 

(3) The history of gas within dark-matter haloes. When a 
halo first forms, it is assumed to have au isothermal- 
sphere density profile. A fraction f2(,/0 of this is in the 
form of gas at the virial temperature, which can cool 
to form stars within a single galaxy at the centre of 
the halo. Application of the standard radiative cooling 
curve shows the rate at which this hot gas cools below 
10* K, and is able to form stars. Energy output from 
supernovae reheats some of the cooled gas back to the 
hot phase. When haloes merge, all hot gas is stripped 
and ends up in the new halo. Thus, each halo main- 
tains an internal account of the amounts of gas being 
transferred between the two phases, and consumed by 
the formation of stars. 

(4) Quiescent star formation. This is one of the more dif- 
ficult areas to model. Most authors assume a star- 
formation timescale that is linearly proportional to the 
circular velocity of the dark-matter halo in which the 
galaxy sits. The star formation rate is equal the ratio of 
the amount of cold gas available and the star-formation 
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timescale. The amount of cold gas available depends 
on the merger history of the halo, the star formation 
history, and the how much cold gas has been reheated 
by feedback processes (discussed below). 

(5) Starbursts. We also introduce the idea that the star- 
formation rate may suffer a sharp spike following a ma- 
jor merger event. This is motivated empirically by the 
existence of ultra-luminous IRAS galaxies, and it allows 
hierarchical models to yield behaviour resembling tradi- 
tional monolithic collapse. Our specific starburst model 
is described in Appendix B. 

(6) Feedback from star-formation. The energy released 
from young stars heats cold gas in proportion to the 
amount of star-formation, returning it to the reservoir 
of hot gas. This ingredient is essential to keep low-mass 
galaxies and proto-galaxies from using up all their gas 
at high redshift. 

(7) Stellar evolution and populations. Having formed stars, 
we wish to predict the appearance of the galaxy that 
results. For this, it is necessary to assume an IMF, and 
to have a spectral synthesis code. Our work assumes 
the spectral models of Jimenez et al. (1998); for solar 
metallicity, the results arc not greatly different from 
those of other workers. The IMF is generally taken to 
be Salpeter, but amy choice is possible. Unlike other 
workers, we taJie it as established that the population 
of brown dwarfs makes a negligible contribution to the 
total stellar mass density, and wc do not allow an ad- 
justable M/L ratio, T, for the stellar population. 

(8) Chemical evolution. It would of course be unrealistic 
to assume that the entire star-formation history of a 
galaxy can unfold at constant metallicity. The evolution 
of the metals must be followed, for two reasons: (i) the 
cooling of the hot gas depends on metal content; (ii) for 
a given age, a stellar population of high metallicity will 
be much rodder. The models of Jimenez et al. (1998) 
allow synthetic stellar populations of any metallicity to 
be constructed. Appendix C discusses how the chemical 
evolution of the gas is followed. 

2.2 Model parameters 

This framework is rather general, and includes a number of 
parameters. Note that there are also parameters involved 
in the cosmological model. In this paper, we consider only 
a single variant of CDM {Q = 1, h = 0.5, n = 1, fift = 
0.06, (Tg = 0.62), in order to illustrate the effects of changes 
in the galaxy-formation model. This cosmological model is 
identical to the one used by Cole et al. (1994). 

The equations containing the specific parameterizations 
of the above assumptions are given in Appedix A. The values 
assumed for these parameters in the models investigated in 
this paper are given in Table 1. 

3 THE TULLY-FISHER / LUMINOSITY 
FUNCTION DISCREPANCY 

Matching both the TuUy-Fisher relation, which relates the 
infrared rest-frame magnitude to the circular velocity Vc of 
the halo, and the B-band luminosity function has been a 
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Figure 1. The Tuily-Fishcr relation: observational data and a 
simple model: single redshift collapse, with all available gas form- 
ing stars in an 0.1 Gyr bm'st (solid lines), or eontiimously with 
an exponential deeay-time of 10 Gyr (dotted lines). The squares 
are individual galaxies as observed by TuUy et al. (1998), whereas 
the thick straight line is an average over four fits to observational 
data by four different authors (see text for details). The thick 
curved line is a fit to the model of Cole et al. (1994, curve taken 
from Heyl et al. 1995). Note that the models of CAFNZ do have 
significant star formation up to the present epoch. 

problem for most published phenomenological galaxy forma- 
tion scenarios. In the following we use, as a starting point, 
the models of Cole ct al. (1994, CAFNZ from here on), which 
show this discrepancy as well. When CAFNZ set their free 
parameters to match the B-band luminosity function, they 
did not find a good match to the I-band TuUy-Fisher re- 
lation: it is either a factor of two too large in Vc, or 2-3 
magnitudes too faint in I. 

We believe that this is largely due to an excess of 
bright galax:ies with rather large circular velocities, caused 
by overmerging within the (extended) Press-Schechter for- 
malism used by CAFNZ. Because of the overabundance of 
massive galaxies, CAFNZ were forced to put in place sev- 
eral 'brakes' in order to prevent the luminosity evolution 
from overshooting. Firstly, they adopt a value for the star- 
formation timescale which is half that of the value of 4 
Gyr that is found in the numerical simulations of Navarro & 
White (1993), on which CAFNZ based their star formation 
models, and which was used to set many other parameters. 
Secondly, they adopted a very large mass-to-light ratio for 
the stellar population: they set T = 2.7, where T is defined 
as the ratio of all stars to that of stars above 0.1 Mq. Fi- 
nally, they assume that all hot gas has primordial metallicity 
(even gas reheated by solar metallicity supernovae), which 
means it cools less rapidly than enriched gas and results in 
loss stars being formed per unit mass. They also assumed 
instant enrichment to solar metallicity for their stellar pop- 
ulations, which helps keep the B-band magnitudes faint. 
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Table 1. Parameters for the phenomenological galaxy formation models discussed in Section 2, 3, and Appendix A. 
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140 


100 


feedback power-law index 
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Stellar populations: 
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3.1 Building blocks for the TuUy-Fisher 
relation 

In order to provide a frame of reference for interpreting 

and understanding the TuUy-Fisher relation produced by 
our models, we first greatly simplify the modelling by only 
incorporating some of the ingredients. The very simplest 
model is one in which a galaxy turns all its gas into stars 
at a single epoch of formation, either in a short burst of 
star formation, or as an exponentionally decaying, but con- 
tinuous star formation mode. We will refer to this type of 
model as a 'single redshift star formation model'. The circu- 
lar velocity is calculated using the spherical collapse model 
(e.g. Peebles 1980), where the density at first collapse is 
178 times the background value, which gives L oc Mgas oc 

In Fig. 1 we plot the resulting Tully-Fisher relation for 

such a model, where the single formation redshift (for all 
galaxies) is annotated for each of the model relations. The 
straight thick solid line is an average of fits to four datasets 
(Pierce & TuUy 1992; Mould et al. 1993 and references 
therein; Giovanelli et al. 1997; Mathewson, Ford & Buch- 
horn 1992), whereas the squares show the data of Young 
(1994), in order to give an indication of the typical spread 
around the observed relation. This average relation corre- 
sponds to L oc Vc'^, so very close to that predicted by the 
spherical collapse model. The most obvious lesson to learn 
from this figure is that we can reproduce the Tully-Fisher 
relation simply by forming all galaxies between a redshift 
of 0.5 and 4, with a peak near z = 1 — 2, and turning all 
available gas into stars shortly after the formation event. 
Galaxies will move straight down in the diagram if only part 
of the gas is used up, so late star formation from a fraction 
of the gas produces the same relation as turning all gas into 
stars at early times. 



3.2 Observed luminosity functions 

In the models of CAFNZ the B and K-band luminosity func- 
tion (LF for short) were used to set some of their free param- 
eters, so that they matched well to the observed luminosity 
functions known at the time: the B-band LF of Loveday et 
al. (1992) and the K-band LF of Mobasher et al. (1986). 
Present day datasets show quite a change in the observed 
LF: Zucca et al. (1997) observe a much steeper faint end for 
the B-band LF, while Gardner et al. (1997) and Glazebrook 
et al. (1995) observe a K-band LF that is fainter at the 
bright end, although there is some disagreement at the faint 
end. It must be gratifying to CAFNZ that their models are 
closer to these new results that to the data of Loveday et al. 
(1992) and Mobasher et al. (1986). The significant discrep- 
ancy with the Tully-Fisher relation remains, however, and 
we attempt to resolve that in the remainder of this Section. 

3.3 The galaxy halo population 

3.3.1 Defining haloes 

In most phenomenological galaxy formation recipes pub- 
lished so far the galaxy halo population and its formation 

and merger history are obtained from the extended Press- 
Schechter (EPS) formalism (see the introduction for refer- 
ences and discussion). It is essential to realize that the EPS 
formalism deals only with systems of density approximately 
200 times the mean, which are assumed to be virialized. 
Such sets of particles arc readily identified in N-body sim- 
ulations e.g. by linking particles via the percolation (a.k.a. 
'friends-of-friends') algorithm. Indeed, the term 'halo' is not 
infrequently taken implicitly to denote only systems of this 
sort. However, we feel this is an unfortunate useage, and we 
will use the term in this paper in a more general sense: as a 
virialized local density maximum. A rich cluster is an EPS 
halo, but the galaxies it contains are defined by dark-matter 
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haloes that constitute substructure in the cluster dark mat- 
ter. These systems collapsed before the cluster, and so have 
present-day density contrasts well above 200. As the clus- 
ter forms, the outer parts of these galaxy-scale haloes merge 
and lose their identity; the high-density cores of the haloes 
will nevertheless survive and remain virialized, and mark 
the locations of galaxies within the cluster. The aim of a 
galaxy formation model is to predict the stellar content of 
these 'embedded' galaxy-scale haloes. 

3.3.2 Overmerging 

If one takes the galaxy haloes merger history from the EPS 
formalism, all embedding haloes arc retained, but embedded 
haloes are lost; EPS records only when haloes are incorpo- 
rated into a larger system, but subsequent evolution within 
the new system is not followed. We adopt the term 'over- 
merging' to describe this loss of information concerning the 
substructure in embedding haloes. 

The use of the EPS formalism is usually motivated 
by the fact that the distribution of masses of the haloes 
and their formation histories are in good agreement with 
those found in N-body simulations (Efstathiou et al. 1988; 
Somerville et al. 1998 and references therein). However, this 
result depends on the choice of the group finder adopted to 
identify haloes in the simulations (Suginohara & Suto 1992). 
The usual choice is percolation with a fixed linking-length, 
which shares with the EPS formalism the property that em- 
bedded haloes arc not identified. It also identifies any over- 
dense group, not just virialized haloes. 

Furthermore, overmerging also occurs in the N-body 
simulations to which the EPS formalism is compared. Small 
groups of particles that represent galaxy haloes get dis- 
rupted by numerical two-body heating, especially in large 
overdense systems (Carlberg 1994; van Kampen 1995). If 
the numerical resolution is good enough, this problem is re- 
solved (Moore et al. 1998; Tormen, Diaferio & Syer 1998), 
but both the spatial and mass resolution need to be rather 
high: Klypin, Gottlober & Kravtsov (1997) find that haloes 
survive for resolutions of a few kpc and 10^-10^ Mq respec- 
tively. Unfortunately, for N-body simulations on a cosmo- 
logical scale, this requires the use of far too many particles 
(on the order of 10^) to be practical, so a special N-body 
technique is needed instead. 

3.3.3 Resolving the overmerging problem 

In order to prevent overmerging, we adopt a simulation tech- 
nique in which a galaxy halo formation recipe is added to an 
otherwise standard N-body technique (van Kampen 1995, 
1997). The idea is that a group of particles that has col- 
lapsed into a virialised system with a mass appropriate of 
a galactic halo is replaced by a single 'halo particle'. Local 
density percolation, also called adaptive friends-of-friends, 
is adopted for finding the groups (van Kampen 1997). This 
is designed to identify the embedded haloes that the tra- 
ditional percolation group finder links up with their parent 
halo. A Gaussian filter, with a length equal to the average 
nearest neighbour distance for a Poisson particle distribu- 
tion, is applied to calculate a local density uq. This local 

density is then used to calculate a local linking length for the 

— 1/3 

percolation algorithm, scaling as Hq ' . The linking length 



also scales with the particle mass m as m ' , in order to 
account for the larger gravitational pull of the more massive 
halo particles. The basic linking length (for mean density 
and non-halo particles) is half the Poissonian average near- 
est neighbour distance. 

The virial equilibrium criterion is put into a simplified 
form using the half-mass radius Rh, motivated by Spitzer 
(1969) who found that for many equilibrium systems 



A group is considered to be virialized when this criterion is 
satisfied with a margin of 25 per cent. The new halo particle 
is softened according to the size of the group it replaces (thus 
conserving energy), and gets the position and momentum of 
the ccntrc-of-mass of the original group. Because groups of 
particles can contain halo particles, merging is naturally in- 
cluded in this recipe as well. Note that we do not assume 
that each halo particle hosts a single galaxy. At formation, 
each halo particle may contain several pre-existing galaxies, 
and the dynamical friction argument (ingredient 2 of Section 
2.1) has to be invoked to see whether they merge. However, 
our procedure clearly treats a good part of the general merg- 
ing of galaxies directly, which increases the robustness of the 
results. 

If groups arc sufficiently massive compared to the nu- 
merical resolution of the code, they will avoid spurious nu- 
merical disruption in future evolution of the density field. 
There is therefore no need to replace these systems by a sin- 
gle particle. In fact, doing so would pose numerical problems 
as supermassive simulation particles will interact violently 
with 'normal' N-body paxticles. Therefore, an upper limit 
should be set to the mass of a halo particle. 

However, the local density percolation group-finder we 
adopt here has been designed to identify haloes in the field as 
well as in overdense regions. This means that a smooth clus- 
ter containing no subhaloes will be identified as one massive 
halo, while a cluster many subhaloes will not be identified 
as a single halo. The centre of the cluster and the subhaloes 
will all be identified as separate entities. This implies that 
there will automatically be a maximum for the mass of the 
groups found, depending on the environment. This is dis- 
cussed in more detail in van Kampen (1999, in preparation). 

Profcrrably the upper limit to the halo particle-mass 
should equal this mgiximum group-mass, but as it depends 
on the group finder parameters, it is haxd to predict its 
value. The results should be robust with respect to reason- 
able variations in this number, so the actual choice is not too 
important. Wc set the upper mass limit to be 7 x 10"^^7i~"^ 
M0, or about 1600 of the initial simulation particles. 

3.3.4 Our halo population 

We performed two high-resolution simulations, using the 
methods described above, with a (spherical) volume of 
(23.2/i~^Mpc)^, using around 8 x 10'"' particles. The first 
simulation started from initial conditions contrained to form 
a galaxy cluster in the centre of the simulation volume, 
while the second simulation started from unconstrained ini- 
tial conditions. We will refer to the first simulation as the 
'cluster model', and to the second one as the 'field model'. 

The importance of overmerging is most clearly illus- 
trated by looking at cumulative dark halo mass functions 
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Figure 2. A comparison of galaxy halo cumulative mass distributions. Dashed lines show the mass function for the Press-Schechter 

formalism for standard CDM with two different normalizations: erg = 0.67 (thick dashed line), and the COBE normalization (erg = 1.19, 
thin dashed line). The thick dotted line represents the mass function for all haloes identified in a standard dissipationless N-body 
simulation (for trg = 0.67) using the percolation algorithm, while the thin dotted line shows the subset of those haloes that have 
virialized. The thickest solid line shows the mass functions for haloes formed using the recipe of van Kampen (1995, 1997), starting from 
the same initial conditions. Including all haloes found from the percolation algorithm time at the final time results in the mass function 
depicted by the thinnest solid line, while the solid line with intermediate thickness shows the subset of these haloes that have virialized. 
Note the difference between the thin solid line and the thin dotted line, which were obtained using the same group finder. 



as produced by the Press-Schochter formalism, a standard 
N-body simulation, and the present method, as described 
above. 

In Fig. 2 we compare all these to each other. The thick 
dashed line is the halo function from the Press-Schechter 
formalism for the cosmology adopted here, with ag = 0.62. 
The COB£ normalized version (i.e. ag = 1.19) is also shown, 
as a thin dashed line. We ran the field model without the 
galaxy halo formation recipe switched on, i.e. just as a tra- 
ditional N-body simulation, and obtained the mass function 
using the same group finder (see Section 3.3.3) used in the 
full simulation technique. This mass function is plotted in 
Fig. 2 as a thick dotted line. It clearly matches the Press- 
Schechter mass function. However, if wc now plot the mass 
function for virialized haloes only, being haloes which satisfy 
eq. (1) within 25 per cent, it falls short by a factor of two 
(thin dotted line). This is because in the EPS formalism 
all collapsed haloes are assumed to virialize, while haloes 
as identified in numerical sinmlations change through merg- 
ing, secondary infall, tidal forces, etc., and are therefore not 
able to remain virialized at all times. This is even true for 



galaxy cluster haloes (Natarajan, Hjortli & van Kampen 
1997). Thus, virialized haloes form a subset of all haloes at 
any one time. Indeed, the haloes formed in our simulations 
contain about 25 per cent of all the mass in the universe, 
while in the EPS formalism all matter is contained in haloes, 
by design. 

The halo functions that are produced by the full simu- 
lations are shown in Fig. 2 as solid lines. The thickest of the 

three lines represents the halo function for the galaxy haloes 
identified by the recipe. The thinnest solid line shows the 
effect of adding the virialized haloes identified at z = 0, 
including the ones above the upper mass- limit. Note that 
there are only a few more haloes added, so our upper mass- 
limit for halo particles is close to the maximum mass found 
from the halo formation recipe employed. 

The right-most solid line shows the mass function for all 
percolated systems at z = 0: galaxy halo particles, newly 
virialized haloes, and all other collapsed haloes. We see 
that locking matter in haloes at early times has the effect 
that the largest mass haloes that form in both the Press- 
Schechter formalism and standard N-body simulations are 
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now correctly resolved into subhaloes and are not counted 
as single objects. More interestingly, the maximum mass 
found for virialized haloes is close to that found for observed 
galaxies. The core of the simulated cluster is one of these 
most massive virialized haloes. 

Because of our choice for the ujjpc'r-limit to halo parti- 
cles, the most massive haloes axe not treated as single parti- 
cles. However, they are identified in the simulation, and are 
incorporated into the merger tree used for the phenomeno- 
logical galaxy formation model. 

4 RESOLVING THE TULLY-FISHER / 
LUMINOSITY FUNCTION DISCREPANCY 

Before showing our model results, a word of warning: 
throughout this section we plot TuUy-Fisher diagrams that 

contain all model galaxies found from the modelling. We 
should therefore bear in mind that only a subset of all these 
galaxies, i.e. the spirals, actually belong in such a diagram. 
As we do not model the morphologies of the galaxies, we aim 
to at least cover the region of Vc — Mi space occupied by 
the observational data. Galaxies outside that region may in 
fact be regarded as ellipticals, SO's, and irregular galaxies, 
and their properties should be tested using other observa- 
tional relations, like the Faber-Jackson relation, which is the 
analogue of the Tully-Fisher relation for ellipticals. 

One could be tempted to merge the Tully-Fisher and 
Faber-Jackson relations into one diagram, as the circular 
velocity is related to the velocity dispersion through the as- 
sumption of an isothermal density profile. However, the 
Faber-Jackson relation concerns the central velocity disper- 
sion, which is not equal to, and usually larger than, the 
velocity dispersion of the dark matter halo. White (1979) 
found that merger remnants tend towards a density pro- 
file, so if elliptical form through merging, they should have 
declining circular velocity profiles. Observationally, it is very 
difficult to measure halo properties of generally gas-poor el- 
liptical galaxies, whereas spiral galaxies have sufficient cold 
gas at large radii to act as a tracer of the halo potential. 
Thus, we need more detailed modelling before we can at- 
tempt a match to the Faber-Jackson relation, and we there- 
fore restrict ourselves to the Tully-Fisher relation for now. 

4.1 Repeating CAFNZ: model m 

To sec the efi'cct of using a merger tree obtained from numer- 
ical simulations that employ the halo formation recipe of van 
Kampen (1997), instead of a merger tree generated accord- 
ing to the extended Press-Schechter formalism, we applied 
the same prescription of gas dynamics and star formation 
that CAFNZ used, with the same parameters (r^ = 2 Gyr 
etc.). The only other difference with the CAFNZ model be- 
side the halo population is the choice of IMF and stellar 
population models (see Table 1), but we believe these differ- 
ences to be marginal (see Appendix A) . The spherical infall 
model has been adopted to calculate the circular velocity, as 
in CAFNZ. 

The result for the cluster simulation is plotted in Fig. 
3a. What is apparent from this figure is that we reproduce 
the CAFNZ Tully-Fisher relation, but miss the high-Vc end 
of it. Because our model galaxy haloes do not overmerge, the 
large- Vc haloes found by CAFNZ are divided into a number 



of smaller haloes in our modelling (see Section 3.3). The 
effect on the B and K band luminosity funcion is a shift to 
the fainter end, as is shown in Fig. 4a. The thick line is the 
luminosity function found by CAFNZ, while the symbols 
with error bars represent the observational data described 
above. Our model m results quite clearly fit neither the 
observations nor the CAFNZ models. But this is actually 
a very nice result, as it allows us to resolve much of the 
Tully-Fisher / luminosity function discrepancy, as follows. 

4.2 Setting T = 1: model n 

Considering the 'brakes' that CAFNZ had to apply to their 
models, it should come as a relief that our model galaxies 
are too faint in B and K, since the Tband magnitudes were 
already too faint for the Tully-Fisher relation to fit. We 
can thus apply an overall brightening to the modelling in 
order to simultaneously match the Tully-Fisher relation and 
both luminosity functions. This is easily achieved by setting 
T = 1 instead of T = 2.7, which was the value that CAFNZ 
needed to adopt. This brightens the I-band Tully-Fisher re- 
lation by just over a magnitude, with the result that the 
predicted distribution of points now overlaps with the ob- 
servational data (Fig. 3b). If we assume that the brightest 
galaxies in I (for a given Vc) are those that were still forming 
stars up to a fairly recent epoch, and are therefore likely to 
be spirals, the predicted Tully-Fisher relation can be con- 
sidered encouraging, except for the faint end. We discuss 
solutions to that at the end of this Section. 

The luminosity functions are shown in Fig. 4b, in the 
corresponding panel. The luminosity function matches the 
observations very well in the B-band, but falls short in the 
K-band. In the next section we therefore introduce two more 
ingredients that not only render the modelling more realistic, 
but provides a better to the K-band luminosity function: 
chemical evolution and starbursts. 

5 CHEMICAL EVOLUTION AND 

STARBURSTS 

5.1 Adding chemical evolution: model c 

Instead of adopting primordial metallicity for the cooling 
function, and constant solar metallicity for the stellar pop- 
ulations, we now incorporate a chemical evolution model. 
The details of the method we use is described in Appendix 
C. This means that we use a different set of stellar pop- 
ulation models, with evolving metallicities. Once a popu- 
lation forms, the initial metallicity and the star formation 
timescale r* determines what the metallicity will be at later 
times. The enrichment of hot gas by chemically evolving 
stellar populations will increase the number of stars forming 
at later times, because the cooling function depends on the 
metallicity. We assume enrichment to be efficient, i.e. the 
metallicity of the cooling gas is equal to that of the popula- 
tions it cools to (mass-averaged over co-existing populations 
in merged galaxies). 

The effect on the Tully-Fisher relation of the inclusion 
of chemical evolution can be seen in Fig. 3c, whereas the 
effect on the luminosity functions is shown in Fig. 4c. We 
see that the B-band luminosity function ends up too bright, 
whereas the K-band one is still too faint. In other words. 
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Figure 3. I-band TuUy-Fisher relation from our modelling (dots) for various choices of parameters. The first six panels correspond to 
the six models discussed in the text, whereas the bottom left panel shows model b for a subset of the data that should represent an 

observational datasct (sec Section 3.4.6 for details). Wc refer to the Icgcnda and to Table 1 for the exact model specifications. Also 
shown arc fits to the distributions from CAFNZ (curved thin solid line), observational data averaged over four different sets of data 
(straight thick solid line, sec main text for details). In the top left panel wc added observational data by TuUy ct al. (1998), shown as 
squares. Note that the specific redshifts at which the galaxy formation recipe is applied are clearly visible for model b, and to a lesser 
extent for the other models, due to the simplified modelling of starbursts. 
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Figure 4. B and K band luminosity functions from our modelling (thick solid lines), corresponding to the TuUy-Fisher relations 
shown in Fig. 5 (in the same order), and those from CAFNZ (thin solid lines). Please refer to the legenda of Fig. 3 and Table 1 for 
details on the models. Also shown are old observational data (triangles) from Loveday et al. (1992) for the B-band and Mobasher et 
al. (1993) for the K-band, and new observational data (squares) from Zucca et al. (1997) for the B-band and Gardner et al. (1997) and 
Glaaebrook et al. (1995) for the K-band. 



most galaxies end up too blue. The reason for this is that 
chemical enrichment boosts late stars formation with rela- 
tively slowly decaying star formation rates. This means that 
we have to find a way to form stars earlier and with much 
more rapidly decaying star formation rates. This leads us 
to the conclusion that wc need to incoorporate a bursting 
mode of star formation, driven by merging activity, which 
is most significant at early times (e.g. Carlberg 1990). 

Another possibility is the addition of a dust model, 
which will typically render galaxies redder. However, we 
leave this for future work to explore. 



5.2 Models with starbursts only: model s 

In order to sec the effect of adding a bursting mode of star 
formation, we first look at a burst-only model in which con- 
tinuous galaxy formation is effectively switched off. We 
make the burst as strong as possible, with set to a very 
large value, and with a very short decay-time. However, the 
duration of the starburst is set as described in Appendix B. 

The Tully-Fisher relation produced by this pure star- 
burst model is shown Fig. 3d. It is quite similar to the chem- 
ical evolution model, in the sense that more stars are formed, 
which brings enough galaxies on or towards the Tully-Fisher 
relation. The spread is a bit larger, as the earliest galaxies 
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fade more rapidly if they experience just a single starburst, 
therefore ending up fainter than in model c, whereas the 
late mergers are brighter because more gas is used up at the 
present epoch than in model c, which has a much larger star 
formation time-scale. 

However, the luminosity functions, as shown in Fig. 4d, 
are different from the chemical evolution ones in the sense 
that the colours end up much redder, providing a better 
match to the observational data for both bands. 

5.3 Models with both chemical evolution and 
starbursts: models a and b 

The pure starburst model produces fairly monochromatic 
galaxies: all galaxies redden quickly and effectively stop 
forming stars, because most mergers happen at high red- 
shifts (e.g. Carlberg 1990), with a merger rate peaJi at 
a fa 1 (van Kampen 1997). If we now switch on continuous 
galaxy formation again, those galaxies that have some gas 
left (likely to be spirals) , can form a (cosmologically) young 
blue population that broadens the colour distribution and 
moves those galaxies on top of the observed TuUy-Fisher re- 
lation. After all, this relation is based on spiral galaxies in 
which some star formation activity is typically found. In 
order to achieve this we set /b = 100 and r^j, = 0.01 Gyr. 

The resulting Tully-Fisher relation is plotted in Fig. 3e, 
and the luminosity functions in Fig. 4e. They all match 
the observational data reasonably well, except for the faint 
end of the Tully-Fisher relation. However, we can easily 
'lift' this faint end by decreasing the amount of feedback to 
/v = 0.04, which comprises model b. The resulting Tully- 
Fisher relation is plotted in Fig. 3f, and this is clearly shown 
to work. It slightly improves the B-band luminosity function 
at the bright end, but at the same time worsens it at the 
faint end, as shown in Figs. 4f. 

Finally, note that the remaining deviations from the 
observed luminosity functions are likely to be resolved by 
taking into account dust, which will redden the colours of 
the galaxies. 

5.3.1 Observational selection 

So far we have made no attempt at discriminating between 
the model galaxies, which mates quite a difference for the 
Tully-Fisher relation, as only spiral galaxies enter the re- 
lation, and only those for which a reliable circular velocity 
can be measured. This means that spiral galaxies showing 
any sign of interaction are typically excluded from an ob- 
servational sample. Thus, we should disregard any model 
galaxy that underwent a major merging event within the 
last 1-2 Gyr. Furthermore, galaxies that have not been part 
of a major merger event for a long time, will have a much 
reduced star formation rate, and axe therefore not likely to 
possess an obvious disk component. So we should also dis- 
regard galaxies that have formed at fairly high redshift and 
evolved passively up to the present epoch. 

We select an 'observational Tully-Fisher sample' from 
all model galaxies by selecting only those galaxies that 
satisfy 1 Gyr < tj^gt < 6 Gyr, where tj^gt is the look- 
back time for the last major merging (or formation) event. 
This corresponds to (approximately) a redshift range of 
0.05 < «iast < 0.5. These are plotted in Fig. 3g. Clearly, 



this subset of the model data provides a reasonable match to 
the observational Tully-Fisher data. However, this attempt 
at observational selection is rather crude, and will need to 
be refined in future work. A proper disk model is a neces- 
sity, as is a dust model which takes the inclination of the 
disk into account. 

6 DISCUSSION AND CONCLUSIONS 

We have argued that overmerging is an important prob- 
lem for phenomenological galaxy formation models, and is 
largely responsible for the discrepancy found in earlier work 
between the Tully-Fisher relation and the luminosity func- 
tion. We resolved the overmerging problem by using an N- 
body simulation technique that includes a galaxy halo for- 
mation recipe. We combined this technique with simplified 
gas dynamics and star formation physics, and a recent stellar 
population model, in order to phenomenologically describe 
galaxy formation and evolution. We included chemical evo- 
lution, and two modes of star formation: major bursts of 
star formation, and a quiescent mode in which stars form 
more continuously. 

With this set-up we match both the B and K band 
luminosity function and the I-band Tully-Fisher relation, 
for an n = 1 standard CDM structure formation scenario. 
Resolving the overmerging problem is the major contributor 
to this result, but the inclusion of chemical evolution and 
starbursts are also important ingredients. 

The new ingredients we have added to the modelling 
of galaxy formation are needed in order to make the mod- 
els more realistic, and are not introduced simply in order 
to give yet more free parameters. Nevertheless, our resolu- 
tion to the Tully-Fisher / luminosity function discrepancy 
may well not be unique, and various other changes to the in- 
gredients of the phenomenological galaxy formation recipe 
might produce similar results. For example, we have not 
studied the influence cosmological parameters have on the 
model galaxy populations, where f2. A, and as are likely to 
be the important parameters. Other types of ingredients 
are possible as well: Somerville & Primack (1998) resolve 
some of the discrepancy using a dust extinction model plus 
a halo-disk approach to feedback. 

We intend to explore these issues in future work. One 
way of resolving the worries about degeneracies in the cos- 
mological/physical parameter space will be to include data 
at intermediate and high redshifts. Nevertheless, although 
the present work has concentrated on modelling the prop- 
erties of the low-redshift universe only, we believe that the 
issues we have raised are sufficiently general that they will 
invariably be part of any successful model for galaxy forma- 
tion. 
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APPENDIX A: DETAILS OF THE GALAXY 
FORMATION RECIPE 

This Appendix provides some further details of the phe- 
nomenological modelling of galaxy formation adopted for 
this paper. We chose to use the modelling of CAFNZ as 
a starting point, but note that many authors use similar 
models with differerences in the details or in the choice of 
parameters. Significant changes to the CAFNZ models are: 
the halo merger history and mass function, the inclusion of 
star bursts, and the modelling of chemical evolution. All pa- 
rameters encountered in the various ingredients of the model 
are listed in Table 1, along with their values for the specific 
models discussed in the main paper. 

HALO POPULATION AND ITS MERGER 
HISTORY 

In most phenomenological galaxy formation recipes pub- 
lished so far the galaxy halo population and its formation 
and merger history are obtained from the extended Press- 
Schechter (EPS) model (see the introduction for references). 
Instead, we use the haloes found an N-body simulations 
which includes a recipe for the formation of galaxy haloes 
(van Kampen 1995, 1997), as described in Section 3.3.3. 
The actual N-body code used is the Barnes & Hut (1986) 
treecode. 

Haloes are assumed to form according to the spherical 
collapse model (e.g. Peebles 1980), and are assumed to settle 
as isothermal spheres, with a constant circular velocity Vc- 
The spherical collapse model provides a relation between 
the mass of the halo, its circular velocity, and its formation 
redshift (e.g. White 1996) : 



Vc = 



Ml 



halo 



,2.35 X 105/i-iM( 
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1 + Zfo 



kms 



(Al) 



where a truncated isothermal density profile has been as- 
sumed. One is forced to adopt this crude model as the ex- 
tended Press-Schechter formalism that is typically used does 
not provide information on the kinematics of the haloes. The 
N-body simulations do provide us with that information, so 
we can obtain Vc directly from the velocity dispersion of 
the halo, which is in virial equilibrium by definition. Again 
assuming an isothermal density profile, we have 
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Figure Al. Comparison of circular velocity from the spherical 
infall model and from N-body simulations (top panel), and a 
similar comparison for the half-mass radii (bottom panel). The 
solid lines simply indicate equality. 



where {v ) is obtained from the numerical simulation. Ba- 
sically, we allow dark haloes with the same mass and «form 
to have a range of values for R^, and therefore different Vc, 
while CAFNZ assume the fixed value 
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In order to access the significance of the difference be- 
tween using lc({«^)) and Vc(M, z), we plotted Vc{M, z) ver- 
sus Vc({w^)) (Fig. Al). Clearly, on average the spherical 
infall model is a good approximation. The scatter repre- 
sents the spread in half-mass radii and velocity dispersions 
for a given mass found in the simulations. An upper limit 
to the mass of a galaxy halo particle needs to be set, as 
discussed in detail in Section 3.3.4. 

In practice, wc sot the upper mass limit to be 7 x 
10^^ M.Q, or about 1600 of the initial simulation par- 
ticles. Although justifyable in numerical terms, this mass 
is also close to the one that results from the cooling crite- 
rion tcool < icoUi where t^ow is the halo collapse time for 
the spherical collapse model (Rees & Ostriker 1977; Binney 



1977; Silk 1977). For primordial abundances, this limit is 
about 1.5 X lo"nb/»"^M0 (White 1996), i.e. « 1.0h~^ x 
10^^ M0 for our choice of fl^,. For enriched gas, with abun- 
dances near solar, the limit goes up to 5 — 10 x 10^^h~^MQ. 

FORMATION AND MERGING OF GALAXIES 
WITHIN DARK HALOES 

For the modelling of the galaxies that populate the dark 
matter haloes we follow the approach of CAFNZ. In their 
scheme, galaxies form iu the centre of dark haloes. When 
dark haloes merge, their galaxies merge with a delay giving 
by the dynamical friction timescale. During this dynamical 
friction phase a dark halo can thus contain several galaxies 
at the same time. One of these galaxies will be designated as 
the galaxy to which halo hot gas can cool, while the others 
are considered satelites. It is important to realize that this 
part of the CAFNZ prescription is kept identical, whereas 
the formation and evolution of the dark halo population is 
different. Thus, dynamical friction is assumed to delay the 
merging of galaxies within merged haloes, with the delay 
given by 



mrg 



(A/halo/A^gal) 



(A4) 



In the following, all gas dynamics and star formation hap- 
pens in between 'major events', which are either the for- 
mation of a halo (plus galaxy), the merger of one or more 
haloes, or the merger of one or more galaxies within an al- 
ready merged halo. These three events define lifetimes for 
haloes and galEixies over which gas can cool, stars can form, 
etc., according to processes we describe next. 

STARS AND GAS 

Gas dynamics and star formation are modelled using an- 
alytical relations that are power-law functions of the halo 
circular velocity, Vc, which is a constant for a given halo as 
we assume an isothermal density profile. 

Each galaxy has a reservoir of stars, cold gas, and hot 
gas. A newly formed galgixy has all its baryons, with a total 
mass of OfcMiiaioj in the form of hot gas. The gas tempera- 
ture is assumed to quickly settle to the virial temperature, 
which is a function of Vc only: 
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The three constants /i, mp and k are the mean molecu- 
lar weight, proton mass, and Boltzmann's constant, respec- 
tively. Hot gas will cool radiatively at a rate which depends 
on the density p, temperature T, and metallicity Z of the 
hot gas. The cooling time-scale is given by 

_ 3 Pgas(r) fcTgas . 

2 n2(r)A(Tgas,Zgas)' ^ ' 

where rie is the electron density, and A(T, Z) the cooling 
function. For the latter we take the cooling functions for 
primordial and solar metallicities as given by Sutherland 
& Dopita (1993), and obtain the cooling function for any 
given metallicity by linear interpolation or extrapolation of 
log A(T, 0.0002) and logA(r,Z0). 

For gas at a constant temperature with a homogeneous 
chemical abundance, the cooling rate depends only on the 
density, and is therefore a function of radius. Thus, at a 
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large enough radius rcooi the cooling time will be larger than 
the time passed between two major halo events (formation 
or merger) . The mass encompassed by this cooling radius is 
transferred from the hot to the cold gas reservoir between 
these two events. 

Stars can form from cold gas, thus transferring mass 
from the cold gas reservoir to the stellar population of the 
galaxy. In the present modelling, mass ejected by dying stars 
is not transferred back to any of the gas reservoirs. However, 
cold gas can be reheated and thus transferred back to the 
hot gas reservoir by the energy output of supernovae. This 
process is called 'feedback' for short. 

The star formation rate is proportional to the amount 
of cold gas actually available: 

dM,{t, Vc)/dt = Meoid(t, Vc)/t.{Vc) , (A7) 

with the star formation time-scale given by 

MVc) = r°(l/c/300 km s'^f' . (A8) 

It is assumed that the reheating of cold gas by supernovae 
can be modelled as 

dMhot(t, Vc)/dt = /3(Vc)dM*(t, Vc)/dt , (A9) 
where the feedback proportionality parameter (3 is given by 

= (Fc/yhot)""""* ■ (^10) 

The overall feedback parameter fv , defined as the fraction of 
energy released by supernovae that is dumped into the cold 
gas reservoir as kinetic energy, determines the values of the 
constants a*, Qhoti and Vhot- In this paper we mostly use 
fv = 0.2, for which CAFNZ found the following values from 
fits to simulation data by Navarro & White (1993): a* = 
— 1.5, Qhot = 5.5, and Vhot = 140 km s~^. For a detailed 
description of the feedback process and its parameters we 
refer to CAFNZ. 

The merging of haloes and galaxies also determine 
where stars, cold, and hot gas end up. A merger of haloes 
brings together all amounts of hot gas of the merging haloes, 
and reheat the new single hot gas reservoir to the virial tem- 
perature of the newly formed halo. It also brings together 
all galaxies that were part of the original haloes. The most 
massive of these will now be the one to which hot gas can 
cool and form staxs. CAFNZ set the circular velocity of this 
halo to be that of the newly formed dark halo, which can 
be rather large, unless its merger time Tmrg, as given by eq. 
(A4), is larger than the time available until the next merger 
event (or the present epoch), in which case the original Vc 
is retained. 

All other galax;ies retain their identity (their stars and 
cold gas reservoir) if they do not merge with the most mas- 
sive one, i.e. if Tmrg is smaller than the time available until 
the next merger event (or the present epoch). In that case 
stars keep forming from the cold gas left in each of those 
galaxies, whereas feedback not only reheats part of the gas, 
but also expells it to the hot gas reservoir of the common 
dark halo. 

STELLAR POPULATION MODELS 

To predict the spectra and photometric properties of galax- 
ies we have used the extensive library of synthetic stellar 
population models computed by Jimenez et al. (1998; 1999), 



where a detailed description of the models can be found. 
Here we briefly review the main ingredients of the models. 
We computed simple synthetic stellar populations (SSPs), 
i.e. stellar populations with fixed metallicity and formed in 
a burst of infinitesimal duration, for metallicities between 
1/100 Zq and 2 Z© and ages between 0.01 and 14 Gyr. 
The library uses new interior models (Jimenez et al. 1998) 
and a new set of stellar photospheres for the coolest mod- 
els (Tgfj < 6000), while Kurucz (1992) models are used for 
Tgjf > 6000K. Two important new features of the models 
are the inclusion of an accurate treatment of all post-main 
sequence evolutionary stages and the use of better theoreti- 
cal photospheric models for low temperatures. 

Simple stellar populations (SSP's) are the building 
blocks of any arbitrarily complicated population since the 
latter can be computed as a sum of SSP's, once the star 
formation rate is provided. In other words, the luminosity 
of a stellar population of age to (since the beginning of star 
formation) can be written as: 

Lx{to)= [" [ ' Lf^{ZM-t)dZdt {All) 

Jo JZi 

where the luminosity of the SSP is: 

L^^^ {Z,tQ-t) = \ M*{Z,m,t)lx{Z,m,to-t)dm, (AU) 

Jrud 

and l\{Z,rn,tQ ~ t) is the luminosity of a star of mass m, 
metallicity Z and age Iq — t, Z^ and Zf are the initial and 
final metallicities, and rriu are the smallest and largest 
stellar mass in the population and dM^:(Z,m,i) /dt is the 
star formation rate at the time t when the SSP is formed. 

We assume that the initial mass function (IMF) is uni- 
versal and is given by a Salpeter IMF (a; = 1.35). In order to 
take account for the chemical evolution we have computed 
detailed models that link dM* /dt and Z, i.e. M*(Z, t) since 
we have assumed that the m dependence is constant (a con- 
stant IMF). This is described in Appendix C. 

APPENDIX B: STARBURSTS 

Besides a quiescent mode of star formation, we consider a 
bursting mode of star formation associated with a major 
formation or merging event. This star formation mode can 
be the dominant one, especially at early times when encoun- 
ters and mergers are frequent, and most galaxies form. The 
strength of the starburst will depend on the dumpiness of 
the matter distribution that is induced by an encounter or 
merging event, on the duration of the starburst phase, and 
on the amount of cold gas available to form stars, which 
can be substantially enhanced by strong inflows of gas into 
the central region of the galaxy during the final stages of a 
merger event (Mihos & Hernquist 1996). Here we adopt a 
simplified approach, by assuming that the SFR during the 
starburst phase can be modelled using the same scaling re- 
lations as for the quiescent star formation mode, but with 
different star formation and cooling rates. The feedback 
strength is assumed to be the same for both star formation 
modes. Also, we adopt the same Salpeter IMF for both 
modes, a choice that is supported by recent observational 
data (see Elmegreen 1999). 

The duration of the starburst, t\y, is assumed to be equal 
to the dynamical timescale of the new dark matter halo: 
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tb = Tdyn = -Rh/o'v- Using the spherical infall model (eq. 
Al), the virial theorem (eq. 1), and the isothermal profile 
(eq. A2), we find 

tb « 0.024(1 + «forin)"^''^/i"^Gyr . (51) 

Therefore, a single starburst at low redshift lasts longer than 
one at high redshift, but as the frequency is much lower, 
it will still be at high redshifts that starbursts are more 
important. Note that the actual starburst might last for a 
much shorter time than the dynamical timescale, depending 
on the amount of material available in the bulges of the 
merging galaxies (Mihos & Hernquist 1994), so t\, should be 
taken as a maximum. 

The bursting star formation timescale r*_b is assumed to 
scale with the halo circular velocity Vc in the same manner 
as the continuous star formation timescale r*, i.e. 

r*,h{Vc) = r°b(l^c/300 km s"^)"* (52) 

(analogous to eq. (A8) in Appendix A). 

It is also necessary to introduce a 'burst factor' /b which 
models the enhanced cooling of hot gas and stronger inflow 
of cold gas during the star bursting phase. It is expressed 
as a density enhancement over and above the mean isother- 
mal profile that will form after the bursting event. If fy, is 
larger than unity, cooling will be more efficient, so that star 
formation will be enhanced during the bursting phase. 

Thus, during a starburst of duration ty,, star formation 
is enhanced by having more cold gas available to form stars, 
and by having those stars formed at a higher rate. The 
shorter star formation timescale will result in a different 
stellar population as compared to the continuous mode for 
two reasons: firstly and simply due to the different formation 
timescale, but secondly because the chemical evolution of 
the population will be different. 

APPENDIX C: CHEMICAL EVOLUTION 

Most of the phenomenological models published so far make 
the simplifiying assumption that hot gas cools with pri- 
mordial metallicity, while the stellar populations form and 
evolve at a constant solar metallicity. This 'instant' enrich- 
ment scheme helps to produce galaxies that arc red enough 
in order to match observations, but is only realistic for star 
bursting galaxies, not for galaxies that form stars more con- 
tinuously. Furthermore, the hot gas around metal-ejecting 
galaxies is enriched, which will boost cooling at intermediate 
and low redshifts. We therefore include chemical enrichment 
in our modelling, in the following way. 

The chemical evolution of a stellar population can be di- 
vided into two parts: the chemical enrichment rate at which 
the gas is polluted by dying stars and the formation rate at 
which gas is transformed into stars. In reality those two as- 
pects are related, since the more metals are injected into the 
interstellar medium, the more efficient star formation should 
be, and also vice versa (e.g. low surface brightness galax- 
ies have low metallicity and low star formation rates). In 
this paper, we assume that the star formation rate depends 
solely on the circular velocity of the halo (see Appendix A). 
We then compute the enrichment of the interstellar medium 
using accurate nucleosynthesis prescriptions. We refer to 
Jimenez et al. (1999) for a detailed description of the chem- 
ical evolution equations used and the stellar yields. 
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Figure CI. Metallicity evolution for various star formation time- 
scales r*. 



Metallicity: 0.01 




Figure C2. B-K colour for a stellar population with an initial 
metallicity of 0.01, as a function of age and star formation time- 
scale T*. 

In principle, one should compute the hydrodynamical 

evolution of the baryons for each halo and calculate the sur- 
face density of the baryonic disc in order to determine the 
star formation rate as well as the infall rate and the out- 
fiow from supernovae. Since our simulations do not have a 
hydrodynamical description of the baryonic component, we 
proceed in a different way. We use the chemical evolution 
models of Matteucci & Francois (1989), but with a star for- 
mation rate given by eq. (A7), thus implicitely assuming a 
specific disk model. As our star formation rate is similar 
to that of e.g. Heavens & Jimenez (1999) and Jimenez et 
al. (1998), this disk model will be close to the disk model 
of these authors. With this set-up it is possible to obtain a 
surface that describes r* vs. Z{t), where r* is given in terms 
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of Vc by eq. (A8). We show Z(t) for four different values of 
T* in Fig. CI. For a given t* and initial Z, a population will 
evolve with time as described in Appendix A. An example 
of the B-K colour evolution is shown in Fig. C2. 

APPENDIX D: STAR FORMATION BELOW 
THE N-BODY RESOLUTION LIMIT 

The N-body simulation technique has a lower limit on the 
mass of the haloes identified. In the hierarchical picture 
these must have formed through merging of haloes with 
masses below that limit. We thus need to approximate 
the stellar and gas mass resulting from the evolution along 
the merger tree below the numerical resolution limit. Both 
modes of star formation need to be accomodated, with the 
complication that merging of small-scale structure is rel- 
atively frequent at early times. Although both modes will 
occur intermittently, we approximate the average of all these 
periodic contributions taken together, by assuming a basic 
star formation timescale of 0.5 Gyr. We estimate the amount 
of gas cooled in the sub-tree by considering a single popu- 
lation which has two thirds of the circular velocity of the 
halo in which the sub-tree terminates, and half its age at 
the time of formation. 

This approximation will be improved upon in future 
work by linking analytical merger subtrees, obtained using 
the extended Prcss-Schechter formalism, to the loaves of the 
merger tree extracted from the numerical simulation. For 
the present paper, these subtrees below the N-body resolu- 
tion limit are not important, as they represent old stellar 
populations that have faded over many Gyr, and are likely 
to end up in bulges and centres of elliptical galaxies. They 
therefore contribute little to the present day luminosities 
that enter the TuUy-Fisher relation and luminosity func- 
tions. 



